Spatiotemporal patterns in a DC semiconductor-gas-discharge system: 
Stability analysis and full numerical solutions 
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A system very similar to a dielectric barrier discharge, but with a simple stationary DC volt- 
age, can be realized by sandwiching a gas discharge and a high-ohmic semiconductor layer between 
two planar electrodes. In experiments this system forms spatiotemporal and temporal patterns 
spontaneously, quite similarly to e.g., Rayleigh-Benard convection. Here it is modeled with a sim- 
ple discharge model with space charge effects, and the semiconductor is approximated as a linear 
conductor. In previous work, this model has reproduced the phase transition from homogeneous sta- 
tionary to homogeneous oscillating states semiquantitatively. In the present work, the formation of 
spatial patterns is investigated through linear stability analysis and through numerical simulations 
of the initial value problem; the methods agree well. They show the onset of spatiotemporal pat- 
terns for high semiconductor resistance. The parameter dependence of temporal or spatiotemporal 
pattern formation is discussed in detail. 

PACS numbers: 05.45.-a, 52.80.-s, 47.54.-r, 02.60.Cb 
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I. INTRODUCTION 



A. Experiments and observations 

Spontaneous pattern formation is a general feature in 
the natural and technical sciences in systems far from 
equilibrium 1]. It is a fascinating phenomenon, but can 
also be detrimental when homogeneity and stationarity 
are required in a technical process. Pattern formation oc- 
curs frequently in gas discharges, like in dielectric barrier 
discharges @, H that are used, e.g., for ozone produc- 
tion and in plasma display panels. It is therefore both 
fundamentally interesting and technically relevant to un- 
derstand the mechanisms and conditions of spontaneous 
symmetry breaking in such systems. 

A dielectric barrier discharge system consists of a lay- 
ered structure that in the simplest case consists of a pla- 
nar electrode, a dielectric layer, a gas discharge layer, and 
another planar electrode 03- To the outer electrodes, 
an AC voltage is applied that forces the system period- 
ically in time. In the present paper, an even simpler 
physical system will be analyzed: a system with essen- 
tially the same set-up, but with a DC voltage supply, i.e., 
a stationary drive. The system is illustrated in Fig. Q] 
To allow current to flow with the DC drive, the dielectric 
layer is replaced by a high-Ohmic semiconductor. 

Independently of the above technical applications, 
both the AC and the DC system have attracted much at- 
tention in the pattern formation community in the past 
two decades, since they are easy to operate, have con- 
venient length and time scales and a wealth of spon- 
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FIG. 1: Scheme of the layers of semiconductor and gas dis- 
charge sandwiched between electrodes with DC voltage. 



taneously created spatio-temporal patterns. When the 
transversal extension of the layers is large enough, ex- 
periments show homogeneous stationary and oscillating 
modes, and patterns with spatial and spatio-temporal 
structure like s trip es, spots, and spirals 0, H, 0, 0, H, H, 
ES EH El G3 EI EE El El El, El- The aspect ratio 
of a thin layer with wide lateral extension and the ob- 
servation from above are reminiscent of Rayleigh-Benard 
convection as the classical pattern forming system in hy- 
drodynamics. 

The experiments on DC driven "barrier" discharges 
111 Ref. g 0, i, 1, El describe not only phenomena at 
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one particular set of parameters, but in @, [n| also ex- 
plore parameter space and draw phase diagrams for the 
transition between different states; therefore we concen- 
trate on those — we are not aware of other experimen- 
tal investigations of such phase diagrams. The experi- 
ments 0, [HI are performed on a nitrogen discharge at 
40 mbar in a gap of 1 mm. The semiconductor layer 
consists of 1.5 mm of GaAs. The applied voltages are 
in the range of 580 V to 740 V. Through photosensitive 
doping, the conductivity of the semiconductor layer can 
be increased by an order of magnitude; the dark con- 
ductivity is a s = 3.2 x 10 _8 /(Ocm). These parameters 
imply that the product of pressure and distance pd of 
the gas discharge is short, but still sufficiently far on 
the right hand side of the Paschen curve that the tran- 
sition from Townsend to glow discharge is subcritical, 
i.e., that the voltage at stationary operation drops when 
the current rises. The resistance of the semiconductor to- 
gether with the applied voltage constrain the operation to 
this transition regime from Townsend to glow discharge. 
The system forms stationary states of the discharge that 
are homogeneous in the transversal direction, as well as 
spontaneously oscillating states that still are spatially ho- 
mogeneous and also oscillating states that show spatial 
patterns in the transversal directions as well. In math- 
ematical terms, the last two states emerge from the ho- 
mogeneous stationary state through a Hopf-bifurcation 
(leading to spontaneous oscillations) or through a Turing- 
Hopf-bifurcation (leading to spatial and temporal struc- 
tures). While we investigated the Hopf-bifurcation in 
previous work [H|,[23j], we here include possible structures 
in the transversal direction and analyze the occurrence 
of Turing-, Turing-Hopf- and Hopf-bifurcations. 



B. Theoretical approaches and understanding 

Theoretical studies both of the DC and of the AC 
driven barrier discharge are few; they basically fall into 
two classes: (i) simulations or (ii) reduced 2D reaction- 
diffusion models, (i) Simulations can only be carried out 
for particular parameter sets; for these parameters, phys- 
ical mechanisms of pattern formation can be identified 
and visualized, (ii) Reaction-diffusion models 0, [2(| Hl| 
in the transversal 2D plane can be investigated analyti- 
cally, but their derivation from the full 3D physical model 
depends on ad hoc approximations whose range of valid- 
ity beyond the linear regime (20| of the Townsend dis- 
charge remains unclear. 

In previous (22l . l23j and in the present work, we have 
chosen a third line of theoretical investigation. Namely, 
(Hi) we investigate the linear stability of the homoge- 
neous stationary state and identify the physical nature 
of the fastest growing destabilization modes. This allows 
us to derive complete phase transition diagrams. This 
approach is used in many branches of the sciences, but 
has been little explored in gas discharges. The stability 
analysis gives a clue for interesting parameter regimes 



that can be further investigated by simulations. 

In more physical detail, the state of theoretical under- 
standing and simulations is summarized as follows. 

First, for the AC system, the importance of surface 
charges deposited on the discharge boundaries in each 
half-cycle of the voltage drive was identified in [f| and 
elaborated in 0, HH, [2y| . Simulations use the same fluid 
models with self generated electric fields as simulations of 
plasma display panels [13] , most work is performed in two 
spatial dimensions. Fully three dimensional simulations 
have recently been performed in [28|. These simulations 
are in reasonable agreement with experimental results, 
given the uncertainty of the microscopic parameters in 
the discharge models. Note that even in the 2D system, 
simulations do not cover many periods of the AC drive, 
typically not more than 20; therefore, they have to start 
from initial conditions fairly close to the final state and 
cannot follow intermediate transients in detail. 

For the simpler DC driven system, we are not aware 
of studies of the full discharge model coupled to the 
high-Ohmic semiconductor layer, except for a reaction- 
diffusion model of the system in the two transversal di- 
rections, that is developed, e.g., in [l9|,[2(|; as said above, 
this approach does not appropriately resolve the dynam- 
ics in the direction normal to the layers. Of course, 
surface charges on the gas-semiconductor interface again 
have to play a role like in the AC system, as the electric 
field can be discontinuous across the interface [20l. I22I [23j] . 
But a full model needs to account for space charges and 
ion travel times in the bulk of the gas discharge as well, 

cf. m m. 

In our previous work (22l . [23| . we concentrated on 
the purely temporal oscillations that occur in a spa- 
tially homogeneous mode, therefore the analysis was re- 
stricted to the direction normal to the layers, assuming 
homogeneity in the transversal direction. The results 
presented in [22| showed that a simple two-component 
reaction-diffusion approximation for current and volt- 
age in the gas discharge layer is not sufficient to de- 
scribe the oscillations, though such a model is sug- 
gested through phenomenological analogies with pat- 
tern forming systems like Belousov-Zhabotinski reac- 
tion, Raylcigh-Benard convection, patterns in bacterial 
colonies etc. In (23|, we predicted the phase transition 
diagram from a homogeneous stationary to a homoge- 
neous oscillating state. These predictions were in semi- 
quantitative agreement with the experiments described 
in0. 



C. Questions addressed in this paper 

Here we investigate the spontaneous emergence of 
spatio-temporal patterns in DC driven systems, and pre- 
dict in which parameter regimes pattern formation can 
be expected. While our previous work dealt with the sta- 
tionary solutions [30l . |3~D ] or the dynamics in the single 
direction normal to the layers [22I [23| , now the transver- 
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sal direction is included into analysis and simulations as 
well. The experiments described in the paper Q and 
thesis [13] systematically explore the parameter depen- 
dence of pattern formation in the system, and we do 
the same, but theoretically. The system in 0] never 
forms modes that are only spatially structured, but sta- 
tionary, i.e., it never undergoes a pure Turing transi- 
tion. Rather, starting from a homogeneous stationary 
state, either the homogeneous oscillating state from [23[ 
is reached through a Hopf-bifurcation, or a spatially 
structured time-dependent state through a Turing- Hopf- 
bifurcation. Furthermore, for high semiconductor resis- 
tivity, typically a spatially structured oscillating state is 
reached while for low resistivity, the oscillating structures 
are homogeneous in the transversal direction. Within 
the present paper we investigate these transitions first 
through linear stability analysis; furthermore interesting 
parameters regimes are investigated through simulations 
as well. We use a fluid model for gas discharge and semi- 
conductor layer coupled to the electrostatic Poisson equa- 
tion; in the gas discharge model electrons are adiabati- 
cally eliminated to reduce computational costs. 

The paper is organized as follows. In Sec. [TT| we in- 
troduce the model, perform dimensional analysis, and 
reduce the model by adiabatic elimination of electrons. 
In Sec. IIII1 the problem of linear stability analysis of the 
homogeneous stationary state is formulated, the equa- 
tions are rewritten and numerical solution strategies are 
discussed. Sec. IIVI contains the results of the stability 
analysis, first the qualitative behavior of the dispersion 
relation as a function of the transversal wave number fc, 
and then predictions on the parameter dependence of the 
dispersion relations. Sec. |V] presents numerical solutions 
of the full initial value problem and a comparison with 
the stability analysis results; they reveal that both meth- 
ods can be trusted. Finally, Sec. I VII contains discussion 
and conclusions. The numerical details for the solution 
of the full p.d.e. system are given in the appendix. 



II. THE MODEL 

In this section, the simplest model for the full two- 
dimensional glow discharge-semiconductor system is in- 
troduced. Its schematic structure is shown in Fig. 1. For 
the gas discharge, it contains electron and ion drift in the 
electric field, bulk impact ionization and secondary emis- 
sion from the cathode as well as space charge effects. The 
semiconductor is approximated with a constant conduc- 
tivity. The same physics was previously studied, e.g., in 
[12, HH and in our previous papers [22l. \2§i, [30L l3l| . How- 
ever, in these previous papers, any pattern formation in 
the transversal direction was excluded and only the sin- 
gle dimension normal to the layers was resolved. The 
model then only allows for stationary solutions [3(J HH 
or temporal oscillations [12, [H| . We now study the onset 
of patterns in the direction parallel to the layers. If the 
layers are laterally sufficiently extended, there is rotation 



and translation invariance within the plane. Linear per- 
turbations can then be decomposed into Fourier modes. 
These Fourier modes can be studied in a two-dimensional 
setting, i.e., in one direction normal and one direction 
parallel to the layers. They are the subject of the paper. 

A. Gas-discharge and semiconductor layers 

The gas-discharge part of the model consists of conti- 
nuity equations for two charged species, namely, electrons 
and positive ions with particle densities n c and n+: 

dtn c + V • r o = source, (1) 
dtn+ + V • r + = source, (2) 

which are coupled to Poisson's equation for the electric 
field in electrostatic approximation: 

V-E g = — (n+-n ), E g = -V$. (3) 

Here, $ is the electric potential, E s is the electric field in 
the gas discharge, e is the elementary charge, and Eq is the 
dielectric constant. The vector fields T c and T + are the 
particle current densities, that in simplest approximation 
are described by drift only. In general, particle diffusion 
could be included, however, diffusion is not likely to gen- 
erate any new structures, but rather to smoothen out the 
structures found here. The drift velocities are assumed to 
depend linearly on the local electric field with mobilities 

Me > 

T c = -fi c n c ~E g , T + = fi + n + Eg, (4) 
hence the electric current in the discharge is 

J g = e (r+ - T e ) = e (u+n + + fi c n c ^j E g . (5) 

Two types of ionization processes are taken into account: 
the a process of electron impact ionization in the bulk of 
the gas, and the 7 process of electron emission by ion im- 
pact onto the cathode. In a local field approximation, the 
a process determines the source terms in the continuity 
equations §T§ and 

source = \T C \ a (|E ff |) , a (|E S |) = a a • (6) 

We use the classical Townsend approximation 

a(|E|/£ )=exp(-£ /|E|). (7) 

The gas discharge layer has a thickness d g , where sub- 
scripts j or s refer to gas or semiconductor layer, re- 
spectively. The semiconductor layer of thickness d s is 
assumed to have a homogeneous and field-independent 
conductivity a s and dielectric constant e s : 

J s = a s F, s , g = £ s £ V-E s . (8) 

The space charge density q inside the semiconductor with 
constant conductivity is assumed to vanish. 
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B. Boundary conditions 



C. Dimensional analysis 



In dimensional units, X parametrizes the direction par- 
allel to the layers, and Z the direction normal to them. 
The anode of the gas discharge is at Z = 0, the cath- 
ode end of the discharge is at Z — d g , and the semicon- 
ductor extends up to Z = d g + d s . (Below, in dimen- 
sionless units, this corresponds to coordinates (x, z) and 
z = 0, L, L + L s .) 

When diffusion is neglected, the ion current and the 
ion density at the anode vanish. This is described by the 
boundary condition on the anode Z = 0: 



r+ (x,o,t) = o 



n+ (X,0,t) = 0. 



(9) 



The boundary condition at the cathode, Z = d g , de- 
scribes the 7-process of secondary electron emission: 

|r c (x,d g ,t)\= 7 \T + (x,d g ,t)\ 

=>■ (j, e n e (X,d g ,t) = 7M+n+ (X,d g ,t) . (10) 

Across the boundary between gas layer and semiconduc- 
tor layer, the electric potential is continuous while the 
discontinuity of the normal electric field is determined 
by the surface charge 



£ = (esEq'Es - £ E 9 ) • n, 



(11) 



where n is the normal vector on the boundary, directed 
from the gas toward the semiconductor, i.e., in the direc- 
tion of growing Z. The change of surface charge in ev- 
ery point (X, t) of the line Z — d g is determined by the 
electric current densities in the adjacent gas and semi- 
conductor layers as 



as= (j fl - J.) n, 



(12) 



where 3 g and J s are the current densities at Z = d g ± 
in gas and semiconductor. 3 S is given in Eq. ©, J g on 
the boundary due to condition (| 10[) is 



J g Z = 9 (1 +7) e/j, + n + B g . 
Equations (| 1 1 p — ff]~5|) are summarized as 

£ = ( £ S £ E S - £ E„ ) • n 



(13) 



(14) 



= S|t=o + J dt {{I + 7)en+/i+E g - cr s E s ) ■ n. 

This boundary condition is valid in any point (X, t) of 
the gas-semiconductor boundary Z — d g . 

Finally, a DC voltage Ut is applied to the gas- 
semiconductor system determining the electric potential 
on the outer boundaries 

$(X,0,t) = 0, ${X,d g + d s ,t) = -U t . (15) 

Here the first potential vanishes due to gauge freedom. 



The dimensional analysis is performed essentially as in 
[22I [231 [29l [30l . [3TI ] . However, as it is useful to measure 
the time in terms of the ion mobility rather than the 
electron mobility, we introduce the intrinsic parameters 
of the system as 



(m) 



1 



to 
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(1 



1 



r = — , go = £oao^o- (16) 



a 



Here time immediately is measured in units of , while 
in (23|, first the time scale t was used. The intrinsic 
dimensionless parameters of the gas discharge are the 
mobility ratio fi of electrons and ions and the length ratio 
L of discharge gap width and impact ionization length: 



Me r 
The dimensionless time, coordinates and fields are 
R t 



(17) 



(*0' 




o-(r,r) 



e n + (R, t) 



S(r,t) 



e n c (R, t) 
E (R, t) 



(18) 



qo E Q 

Here the dimensional R is expressed by coordinates 
(X, Z) and the dimensionless r by (x,z). 
The total applied voltage is rescaled as 

- (19) 

The dimensionless parameters of the semiconductor are 
conductivity a s and width L s : 



fJ>+Qo 



L s = 



ro 



(20) 



Note that the dimensionless conductivity is now also 
measured on the scale of ion mobility. Dimensionless 
capacitance and resistance of the semiconductor and its 
characteristic relaxation time are expressed in terms of 
these quantities as 



K„ = 



C s = 



= C s Ks = 



(21) 



D. Adiabatic elimination of electrons and final 
formulation of the problem 

The dynamics of a glow discharge takes place through 
ion motion where the ions are much slower than the elec- 
trons. As in [23j], the electrons therefore equilibrate on 
the time scale of ion motion and can hence be adiabat- 
ically eliminated: After substituting s = cr/fx, the gas 
discharge part of the system has the form 

Hd T s-V-(£s) = s\S\a(\S\), (22) 
d T p + V-{£ P ) - s|£|a(|£|), (23) 
V-£ = p-fis, £ = -V0, (24) 
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and in the limit of p — * 0, it becomes 

- V- (£s) = s|£|a(|£|), (25) 
d T p + V • (£p) = s|f|a(|f|), (26) 
V£ = p, £ = -V0. (27) 

As in [23], the electric field £ is now only influenced by 
the ion density p, and not by the much smaller density of 
fast electrons (since the electrons are generated in equal 
numbers, but leave the system much more rapidly), and 
the electrons s follow the ion motion instantaneously: 
given the electron density on the cathode s(x,L,t) and 
the field distribution £ in the gas gap, the electron den- 
sity is determined everywhere through (|25p . The bound- 
ary conditions © and (JTUJ) for electrons and ions are 

pOr,0,T)=0, s{x,L,t)= 1P {x,L,t). (28) 

After rescaling, the semiconductor is written as 

V-£ = 0, £ = -V<(>, j s = <r s £, (29) 

and the condition (fT4")) for the charge on the 
semiconductor-gas boundary becomes 



— i— = [e s £\ z _ L+ -£\,- L - 
q Q /r V ,Z - L+ ,Z - L 



n 



(30) 



qo/r 



r=0 JO 



(It 



((1+7) p£\ 



<j s £\ L+ I ■ u- 



In the perturbation analysis, the differential form of 
charge conservation on the boundary is used 



qo/r 



= ((l + 7 ) p£\ L _ -a s £\ 



L+i 



n. 



(31) 



The total width of the layered structure is L z = L + 
L s . On its outer boundaries z = and z = L z , the 
electrodes are on the electric potential 4> (x, 0, r) = and 
<fi (x, L z ,t) — —Lit, respectively. Finally, in the numerical 
solutions of the PDEs, lateral boundaries at x = and 
X = L x with periodic boundary conditions for <p, p, and 
s are introduced. 



E. Parameter regime of the experiments 

The parameters are taken as in the experiments @, G3 
and in our previous work [23| . The discharge is in nitro- 
gen at 40 mbar in a gap of 1 mm. The semiconductor 
layer consists of 1.5 mm of GaAs with dielectric con- 
stant e s = 13.1. The applied voltages are in the range 
of 580-740 V. Through photosensitive doping, the con- 
ductivity of the semiconductor layer could be increased 
by an order of magnitude; the dark conductivity was 
a s = 3.2 x 10" 8 /(ftcm). 

For the gas discharge, we used the ion mobility 
/(+ = 23.33 cm 2 /(Vs) and electron mobility p e = 
6666.6 cm 2 /(Vs), therefore the mobility ratio is p = 



p+/p e = 0.0035. For a = Ap = [27.8 pm]' 1 and for 
Eo = Bp = 10.3 kV/cm, we used values from [33]. The 
secondary emission coefficient was taken as 7 = 0.08. 
(Comparison with experiment in [23| as reproduced in 
Fig. [3] below might suggest 7 = 0.16, but that seems 
unreasonably large.) Therefore the intrinsic scales from 
(TT^l) are 



ro = 27.8 pm, 
q Q = 2.04- 10 12 e/cm 3 , 



= 11.6 ns, 
E = 10.3 kV/cm, 



(32) 



the gas gap width of d — 1 mm corresponds to L = 36 
in dimensionless units, and the semiconductor width of 
d s = 1.5 mm to a dimensionless value of L s = 54. The 
dimensionless applied voltages are in the range of 17.5 < 
Ut < 50, which correspond to the dimensional range of 
500 V < Ut < 1425 V. The dimensionless capacitance of 
the semiconductor is C s = 0.243. We investigate the con- 
ductivity range of 6- 10~ 8 /(ft cm) < 5 S < 6-10" 7 /(f7 cm) 
which corresponds to a semiconductor resistance 1Z S of 
700 to 7000 in the new dimensionless units (or to 2 • 10 5 
to 2 • 10 6 in the units of our previous papers [13, HH). 

For the lowest conductivity ofa s = 3.2x 10~ 8 /(fi cm), 
pattern formation consistently occurs neither in our anal- 
ysis nor in the experiment; this case is not discussed fur- 
ther in this paper. 



F. Experiment: between Townsend and glow 

The parameter regime of the experiments @, E3] as dis- 
cussed above can now be placed in the transition regime 
between Townsend and glow discharge as follows. We 
recall [30l [3l[ that the stationary voltage Wxown of the 
space charge free Townsend regime is minimal at a dis- 
charge length of L = e In [(1 + , which is L = 7.07 
for 7 = 0.08 or for example L — 18.8 for 7 = 10~ 3 . 
(^Town(^) is known as the Paschen curve.) Continu- 
ing with 7 = 0.08 in the remainder of the paper, the 
transition from Townsend discharge to the space charge 
dominated glow discharge is purely subcritical, i.e., the 
current-voltage characteristic is falling, if the discharge 
is longer than L cr n = e 2 In [(1 + 7V7] = 19.23; this is 
here the case. In fact, for the experimental value L = 36 
and 7 = 0.08, the Townsend voltage is Wtowd = 13-7, 
while the minimum of the voltage in the glow regime is 
Ug\ OV j = 11.5, it is reached at a current of J w 0.3 as 
shown in figure [21 In the experiment the applied volt- 
age does not exceed bit — 50, and the resistance of the 
semiconductor is at least 1Z S — 700, this situation is in- 
dicated as the dashed load line hi = hi t — 1Z S J in fig- 
ure Therefore the current in the stationary homoge- 
neous mode does not exceed 0.07, i.e., it stays in the 
transition region between Townsend and glow regime. 
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FIG. 2: Solid line: Current voltage characteristics of the gas 
discharge for L — 36 and 7 = 0.08. The Townsend voltage 
Wtowii = 13.7 has vanishing current J = 0, the minimal volt- 
age in the glow regime is W g i ow = 11.5 at J w 0.3. Dashed 
line: load line of the semiconductor for Ut = 50 and 1Z S = 700. 
For smaller Ut and larger TZ S , the load line moves closer to 
the U axis, i.e., to the Townsend regime. The intersection of 
(solid) gas characteristics and (dashed) load line indicates the 
stationary solution of the system. 



III. STABILITY ANALYSIS: METHOD 

In this section, the stability analysis of the homoge- 
neous stationary state is set up. While in earlier work 
[HI, only temporal oscillations were admitted, here the 
stability with respect to temporal and spatial perturba- 
tions is analyzed. In particular, linearized equations are 
derived that define an eigenvalue problem, and the nu- 
merical solution strategy is discussed. It becomes partic- 
ularly demanding for large wave numbers k. The method 
developed in this chapter forms the basis for the deriva- 
tion of dispersion relations in Chapter HVl 



Linear perturbation analysis for transversal 
Fourier modes: problem definition 



Here the equations are derived that describe linear per- 
turbations of the stationary state that is furthermore 
homogeneous in the transversal direction, in agreement 
with the external boundary conditions. The analysis is 
set up as in [23|, but now also transversal perturbations 
are admitted. The unperturbed equations are denoted 
by a subscript 0, they are for p — > 



d z jo = -joa(£o), 
£oPo+jo = Jo, d z J 
d z £o = po, d z cj) 



where j = s Q £ , (33) 
= 0, (34) 
= -So, (35) 



with boundary conditions 

io(0) = Jo, Mo) ■ 
jo(L) = — ^ J , (j) (L) 
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(36) 

-W , Ut=U +K s J . 



Here Jo is the total current and U t is the applied voltage. 
For a further discussion, see [H, IS Hl| . 

The first order perturbation theory is denoted by a sub- 
script 1. As the transversal modes can be decomposed 
into Fourier modes with wavenumber k within linear per- 
turbation theory, the ansatz 



s(x,z,t) = s (z) + Si (z) e tkx+x \ 
p(x, z, t) = p (z) + Pl (z) e tkx+Xr , 
<I>(x,z,t)=M*)+Mz) e lkx+XT . 



(37) 
(38) 
(39) 



is used where the perturbation is supposed to be small. 
Insertion of this ansatz into the equations for the gas 
discharge lpg |) -([27 |) yields 



d z si = - a(£ ) + 



d z £ 
£0 



so 

si - 7-P1 



d z pi 



d z £i 

d z <j)i 



so lp s 

t-0 



a(£ )si - 

+ (|«(&) 

pi - k 2 4>i. 



d z s 



£0 

A + po + d z £o 



s a'(£ ) ) £1 



£0 
dzPo 
£0 



Pi 



+ s a'(£ ) £j_ 



-£1. 



(40) 



(41) 



(42) 
(43) 



Here £ 1 is the field perturbation in the z direction. The 
boundary conditions are: 



pi(0)=0, <M0)=0, S1 {L)= 1P1 {L), 



(44) 



where z = L is the boundary between gas discharge and 
semiconductor layer. 

In the semiconductor layer, the equation = with 
the boundary condition <f>(L z ) = —hit at the position of 
the cathode L z = L + L s has to be solved. For (pi, this 
means that we have to solve A0i = with 4>i(L z ) = 0. 
This problem is solved explicitly for L < z < L z by 



0i(z) = Ci sinh(fc(z - L z )), 



(45) 



with the arbitrary coefficient C\. The 'jump' condition 
(|3"0"|) . (|3T|) for the electric field on the semiconductor gas 
discharge boundary is expressed as 



-C\k cosh(A;J s ) [A e s + cr s ] 

= (1 + 7)(/?o£i + £oPi) + A£i 



(46) 



after E is eliminated. As the potential (|43| is continuous 
we get on the boundary z — L~ of the gas discharge 
region 



L (L+) = 4> X {L-) = -Ci sinh(A:L s ). 



(47) 
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Now Ci in (gH) can be substituted by fa (L) through ((47 
The result is the second boundary condition at z = L 



fa(L) = 



(1 + j)( Po £i + Sopi) + X£i 



A e s 



tanh(/cL s ) 



(48) 

Now the semiconductor is integrated out, and we are 
left with four first order ordinary differential equations 
(J40])— (J43J) and four boundary conditions (|44]), 63) that 
together determine an eigenvalue problem for A = A(fc). 



B. New fields lead to compacter equations 



It is convenient to write the equations (I40[) in terms of 
the fields h and g 

h = — + -i and g = £ £i (49) 
so So 

as in [23| . Furthermore, for non- vanishing wave-numbers 
k, it is convenient to use charge conservation 



= d T p + V • [s£ + p£j = V • [d T £ + (s + p)£J (50) 

to eliminate particle densities completely in favor of the 
z component of the total current density 

ji = A£i + (si + pi)£ + (s + p )£i. (51) 

This leads to a compacter form of the system (|4"0"]) -([4"3 |) : 



a 1 k 2 

d z h = -^g-'Lfa, (52) 
to to 

d z g = -jo h - — g + ji - k 2 £ fa, (53) 
oo 



d z ji = ~k 2 yX + ^ | <>,. 

dzfa = -7-5- 
co 



(54) 
(55) 



The boundary conditions (|4~4"|) and (|4"B")) are rewritten as 
A 



Ji(0) 



,g(0) + Jo MO), 



£o(0) 

0i (0) = 0, 

h(L) - + Jo h(L), 

1Z S ji(L) tanh(A:L s ) 



ML) = 



1 + A t , 



kLc 



Note that in the last equation, the identity (A e s + <r s ) — 
(1 + \t s )L s /1Z s was used. Note further that the limit of 
k — > of these equations reproduces the limit of /i — > 
of the analogous ID equations in [23|. 



C. Formal solution and numerical implementation 

In matrix form, the linearized equations (|52p - (|55[) arc 



d z v = Ma • v, where v(z) 




(60) 



and M x (z) 



-a'/£ -k 2 /£ a 
-jo -A/f 1 -fc 2 ^o 

-fc 2 (A + J /£ ) 
-l/£ 



The matrix Mj(z) depends on wavenumber k and eigen- 
value A and total current Jo- It depends on z through 
the functions £q(z), a(£o(z)), and jo(z). 

The boundary conditions ([56]) . (|57|) at z = mean that 
the general solution of the linear equation can be written 
as a superposition of two independent solutions of (|60|) 



v(z) = Cl v^(z) + c 2 v^(z), a < vW = M r vW, 
1/Jo \ / 



(0) 



,(2) 



(0) 




(61) 



The solution (j6"Tj) has to obey the boundary conditions 
([58]) and ([59]) at z = L as well. Denoting the components 
of the solutions as vW = f/jW, gW, 0^ ) for i = 1,2, 
we get 



ci 



(i) 



£ - Jo 
to 



-C-2 



(2) 



f <? (2) - Jo 

co 



= 0, 



z=L 



Cl 



-C2 



(1) _ (1 + At.) fcL a (i) 
s Jl tanh(fci s ) 01 

^ sj (2)_(l + Ar^fei £0 ( 2) 



tanh(/cL s ) 



(62) 



0. (63) 



z=L 



These equations have nontrivial solutions, if the determi- 
nant 



AO) 



(64) 



(56) 




3i ] - 




A 2 ) A 12) T 


few 


(57) 






(1+Ar s ) kL s ,(1) 
tanh(fcL s ) V l 


v .(2) (1+Ar s ) fcL 
/V * tanh(feL s ) 




(58) 


vanishes at 


Z = L 






(59) 






A(z = L) 


= 0. 


(65) 



Now for a fixed fc, we start with some initial estimate for 
the eigenvalue A(fc) and solve equation (fBTj) numerically 
for both initial values. The next estimate for A can be 
found from condition (|65p since it is quadratic in A. This 
process is iterated until the accuracy is sufficient. We 
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used the condition |A(" +1) - \W\ / |A( ,l+1 )| < 1(T 8 for 
the n'th iteration step to finish iterations. For the sta- 
bility of the iteration process under-relaxation was used. 
For the integration of the equations (fBTJ)) . we used the 
classic fourth-order Runge-Kutta method on a grid with 
500 nodes. The majority of investigated solutions of the 
present problem are oscillating, therefore the eigenvalues 
A are complex, and the eigenfunctions v(z) are complex 
as well. We have taken this into account by working with 
complex fields. 

After the eigenvalue A(fc) is found, the eigenfunction is 
determined by inserting the ratio 



£2 

Cl 



tanh(fc£ s ) K s j[ 1] (L) - kL s (1 + Ar s ) <f)\ L, {L) 



tanh(fc£ s ) TZ S j{ 2) (L) - kL s (1 + Ar s ) (j>\' J (L) 

(66) 



(2), 



into Eq. 



D. Numerical strategy for large k 

The above method gives reliable results for small val- 
ues of k and has been used to derive the results presented 
in Sect. IIV Bl However, for investigating possible insta- 
bilities for large wave number k as in Sect. IIV A} e.g., for 
finding both solution branches for k > 0(10°) in Fig. 2J 
a better strategy is needed. 

There are two points where the integration routine has 
to be improved for large k. There is first the fact that 
the matrix of coefficients is poorly conditioned. This can 
be seen by noting that one column is much bigger than 
another one. A more precise measure of numerical ill- 
conditioning of a matrix is provided by computing the 
normalized determinant of the matrix. When the nor- 
malized determinant is much smaller than unity, the ma- 
trix is ill-conditioned. The normalized determinant is 
obtained by dividing the value of the determinant of the 
matrix by the product of the norms of the vectors forming 
the rows of the matrix. 

The second point is the so-called 'build-up' error. The 
difficulty arises because the solution (|61[) requires com- 
bining numbers which are large compared to the desired 
solution; that is vW(z) and -v^(z) can be up to 3 or- 
ders of magnitude larger than their linear combination, 
which is the actual solution. Hence significant digits are 
lost through subtraction. This error cannot be avoided 
by a more accurate integration unless all computations 
are carried out with higher precision. Godunov [HJ pro- 
posed a method for avoiding the loss of significance which 
does not require multi-precision arithmetics and which is 
based on keeping the matrix of base solutions orthogonal 
at each step of the integration. 

A modification of Godunov's method [35|, which is 
computationally more efficient and which yields better 
accuracy, is implemented in the algorithm used for large 
k. The main difference to the algorithm described in 
Section fill Cl is that here we examine the base solutions 



(obtained by any standard integration method) at each 
mesh point and when they exceed certain nonorthogo- 
nality criteria we orthonormalize the base solution. Wc 
have to start with initial conditions that are orthogonal 
to each other and not only to the boundary conditions: 



v«(0) 



v< 2 >(0) 



T 





£0(0) 



1 





\ 



where H: 



£0(0) 







(Jo 



£0(0) 



1 + Jo 



£o(0) 



(67) 



Based on the orthogonalization developed in [35| , the dif- 
ferential equation (|60p can be solved very accurately even 
for large k, which allows us to find the eigenvalues accord- 
ing to the criteria described in the previous section. 

Since we must solve the matrix equation (|6"0|) several 
times for different values of A, we must insist that the 
orthonormalization is the same for all these solutions. In 
essence we must insist that the determinant is uniformly 
scaled in order for the successive approximations for A 
to be consistent. Numerically this can be accomplished 
by determining the set of orthonormalization points and 
matrices for the solution corresponding to the first initial 
guess for A and thereafter applying the same matrices at 
the corresponding points for all successive solutions. 

The program is written in C, and the integration 
method is one-step simple Runge-Kutta and on the do- 
main L=36 the number of grid points is varied from 
n=500 to n=18000 depending on the range of k. 



IV. STABILITY ANALYSIS: 
RELATIONS 



DISPERSION 



Having set the mathematical basis for the stability 
analysis, now dispersion relations and bifurcation dia- 
grams can be derived and discussed. In previous work 
[23j, we have analyzed pure Hopf transitions where the 
homogeneous stationary state becomes unstable to ho- 
mogeneous oscillations. The bifurcation diagram for 
the Hopf transition for the experimental parameters de- 
scribed in Sections III El and III Fl is shown in Fig. [31 it is 
reproduced from Fig. 11 in 23]. Any structures in the 
transversal direction were excluded in this earlier work. 
We now take this diagram as a basis and investigate 
which additional spatial or spatio-temporal instabilities 
can occur. 
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FIG. 3: The Hopf bifurcation lines where the homogeneous 
stationary state becomes unstable to homogeneous temporal 
oscillations; while structure formation in the transversal di- 
rection is excluded. The bifurcation is drawn as a function 
of Ut and 1/R S while all other parameters keep the constant 
values described in Sect. Ill El We recall that 1Z S can be varied 
by a factor of 10 by photo- illumination. The oscillations occur 
above the lines; shown are two calculated lines for 7 = 0.08 
and 7 = 0.16, and the experimentally measured line from Q. 
No parameters have been fitted. The figure reproduces Fig. 11 
from 23], therefore the old convention for R s = 'R. S /[J is used: 
The upper axis label 1/R S = 4 • 10 -6 in the plot corresponds 
to the small 1Z S = 875. 



A. Hopf or Turing-Hopf bifurcations 

If patterns form spontaneously in the system due to 
a linear instability, i.e., in a supercritical bifurcation, 
its signature will be found in the dispersion relation 
A = A(fc). More precisely, there will be a band of Fourier 
modes with positive growth rate Re \(k) > 0. If the 
instability is purely growing or shrinking without oscilla- 
tions, i.e., if Im X(k) = 0, it is called a Turing instability. 
On the other hand, if the imaginary part of the disper- 
sion relation does not vanish (Im A(fc) ^ 0), the system 
oscillates. If the most unstable mode k* has no spatial 
structure, k* = 0, but only oscillates, we speak of a Hopf 
transition, while if k* ^ 0, the transition is called Turing- 
Hopf. 

For three values of TZ S , namely for 700, 1400 and 7000, 
the dependence of the dispersion relation A(fc) on the 
applied voltage U t was investigated. 



1. Qualitative behavior for 1Z S = 700 and 1400 

For 1Z S — 700, a generic shape of the dispersion relation 
is presented in Fig. |4] Here a very large range of k values 
is shown on a logarithmic scale. The point k = on the 
axis was previously treated in (23j . The dispersion curves 
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FIG. 4: Real part of the dispersion relation for 1Z S = 700 and 
Ut = 23. 



X(k) extend continuously from k = to small k: there 
is a pair of complex conjugate eigenvalues A(fc) that is 
represented by the single line for Re(A). For k of order 
unity, the two complex conjugate solutions A(fc) in Fig.fi] 
merge and form two solutions with different real part and 
vanishing imaginary part. 

The dynamic behavior is typically dominated by the 
mode with the largest positive growth rate Re(A). If the 
growth rate is negative for all k, the system is dynami- 
cally stable. Here the mode k* with the largest growth 
rate has vanishing wave number k* = and nonvanish- 
ing Im A; therefore we expect a Hopf bifurcation towards 
oscillating homogeneous states. If the upper solution in 
Fig. [4] would develop a positive growth rate for large k, 
we had found an exponentially growing, purely spatial 
Turing instability at short wave lengths, but we haven't 
found such behavior. 

Variation of the applied voltage Ut leads to the same 
qualitative behavior: the temporally oscillating, but spa- 
tially homogeneous mode k = has the largest growth 
rate. Whether this maximal growth rate is positive or 
negative, can be read from Fig. [3] The behavior for 
1Z S = 1400 is qualitatively the same. 



2. Qualitative behavior for TZ a = 7000 

Further increase of the semiconductor resistivity to 
1Z S = 7000 creates a new feature, namely a Turing-Hopf- 
instability: the growth rate becomes maximal for some 
non- vanishing, but very small value of k — k* > 0, as 
can be seen in Fig. [5] 

In Fig. [51 real and imaginary part of the two largest 
eigenvalues are plotted for a larger range of k. Up to k « 
8, a pair of complex conjugate eigenvalues is found, then 
two different branches of purely real eigenvalues emerge, 
similarly to the behavior for smaller 1Z S in Fig. [3J 
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FIG. 6: Real and imaginary part of the two branches of the 
dispersion relation with largest growth rate for the same pa- 
rameters as in Fig. \5\ 



Finally, Fig. [7] shows that the most unstable branch 
of eigenvalues approaches Re A (A;) = from below for 
k —>■ oo, but does not develop any positive growth rate. 
Growth rates for large k always stayed negative when 
exploring the parameters of the system numerically. If 
positive growth rates would exist, they would indicate a 
purely growing spatial mode with very short wave length. 



3. Quantitative predictions 



The above observations are now quantified in Figure [5] 
It shows how the most unstable wave number k* and real 
and imaginary part of the corresponding eigenvalue A(fc*) 
depend on the feeding voltage lit at different dimension- 
less resistances TZ S = 700, 1400 and 7000. The curves 
begin where the applied voltage Ut equals the Townsend 
breakdown voltage 13.7, cf. Section Hi Fl 

Panel (a) shows that for small lit the most unstable 
wavelength k* is always non-vanishing, and that it de- 
creases for growing Lit until it vanishes at some critical 
lit- Panel (b) shows that the growth rate Re A(fc*) of the 
most unstable mode can stay negative in the whole range 
where k* > 0. This explains why the finite wave length 
instabilities were not seen above for small 1Z S . Panel (c) 
shows that the most unstable modes k* are always oscil- 
lating in time. 

These results are further summarized in Figure[9]which 
is the counterpart of Fig. [3l but with 1/TZ S = l/(R s /z), 
i.e., with the definition of 1Z S of the present paper. The 
solid line in Fig. [9] is the line where Re A(fc*) changes 
sign. This line is essentially identical with the solid line 
in Fig. [3] that denotes the sign change of Re A(0). The 
dash-dotted line denotes where the most unstable wave 
number k* vanishes. The dashed lines denote the voltage 
of the Townsend and the glow discharge. 

The figure shows a calculated bifurcation diagram: In 
region I, a discharge cannot form, in region II, the ho- 
mogeneous stationary discharge is stable, in region III, 
it is Turing-Hopf-unstable, and in region IV, it is Hopf- 
unstable. 
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FIG. 8: (a) Most unstable wave number k* , (b) real part and 
(c) imaginary part of the corresponding eigenvalue A(fc*) as 
a function of the total voltage Ut for 1Z S = 700, 1400, 7000. 
The change of 1Z S can be achieved by photo-illumination. 



4- Comparison with experiments 

Comparing this calculated bifurcation diagram with 
the experimental one in 0, fiol ] , there is one point in com- 
mon and one differs. 

The common point is that there are no purely spatial 
patterns without oscillations in the experiment, in agree- 
ment with the theoretical prediction that there is no pure 
Turing instability. Furthermore the transition from sta- 
tionary to oscillating states occurs roughly in the same 
parameter regime, cf. Fig. [3) The oscillations are due to 
the following cycle: (i) the voltage on the gas discharge 
is above the Townsend value and the ionization in the 
discharge increases, (ii) the discharge deposits a surface 
charge on the gas-semiconductor interface hence reducing 
the voltage on the discharge and eventually extinguish- 
ing it, (Hi) the high resistance of the semiconductor leads 
to a long Maxwell relaxation time of the initial voltage 
distribution between gas and semiconductor, after which 
the cycle repeats. 

The difference between experiment and calculation lies 
in the sequence of temporal and spatio-temporal pat- 
terns. In the experiments, for large conductivity 1/TZ S 
on increasing Ut first purely temporal and then spatio- 
temporal patterns are formed. For smaller 1/1Z S , the 




FIG. 9: Calculated bifurcation diagram as a function of volt- 
age Ut and conductivity 1/1Z S - Dashed lines: minimal voltage 
Wgiow = 11-5 (thin) in the glow regime and Townsend voltage 
Wtowii = 13.7 (thick) of the gas discharge according to sec- 
tion III Fl Solid line: Re A(A;*)=0; this line is almost identical 
with the solid line in Fig. [3] that denotes Re A(0)=0. Dash- 
dotted line: k* = 0. Region I cannot form a discharge. In 
region II, the homogeneous stationary discharge is stable, in 
region III, it is Turing-Hopf-unstable, and in region IV, it is 
Hopf-unstable. For future investigations, it is interesting to 
note that all bifurcation lines become almost straight when 
plotted as a function of 1Z S rather than 1/1Z S - 



system goes from stationary to oscillating and back to 
stationary behavior without loosing homogeneity. In our 
calculation, for large 1/TZ S , the system directly transits 
from homogeneous stationary to homogeneous oscillat- 
ing, while for small 1 /1Z S there is first a range of a Turing- 
Hopf-instability, and then a pure Hopf-instability takes 
over. 

On this discrepancy, three remarks are in place. 
(/) The nature of the linear instability of the homo- 
geneous stationary system does not automatically pre- 
dict the fully developed nonlinear pattern. (II) Our 
spatial patterns have wave numbers k* typically much 
smaller than 0.1 which corresponds to wave lengths much 
larger than 60. Therefore our spatio-temporal instabili- 
ties do not correspond to the dynamic filaments described 
in [3], but rather to a range of diffuse moving bands. 
These diffuse moving bands occur before the blinking fil- 
aments discussed in [7j; they are only shortly described 
in the later Ph.D. thesis [10|. We therefore suggest that 
the moving waves or bands can be identified with our 
predicted Turing-Hopf-instability, that creates running 
waves as well; the experimental data in (lfj| do not allow 
a further comparison with theory and we therefore invite 
further experimental investigations. The authors of 

the experimental papers have suggested later [t| Q1J [36[ 
that their experiment is more complex than our simple 
model due to nonlinearities in the semiconductor and in 
gas heating. Our calculation therefore shows that al- 
ready this simple model exhibits Hopf- and Turing- Hopf- 
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FIG. 10: The influence of the width L s of the semiconductor 
layer on real and imaginary part of the dispersion relation 
A(fc), for equidistant L s from the range between L s — 10 and 
150 at Jo = 1.32 • 10~ 5 , L = 36. For L s = 54, all parameters 
are the same as in Figs. [6] and [7] TZ S = 7000, C s = 0.243, 
and U t = 40. 



instabilities in the same parameter range. 



B. Dependence on gap lengths and resistance 

We now investigate how the dispersion relations change 
when other system parameters are varied. While keeping 
material parameters for gas and semiconductor fixed, the 
following physical parameters can be varied: the widths 
L s and L of the semiconductor and the gas layer, the 
externally applied voltage Ut and the semiconductor re- 
sistance 1Z S by a factor of 10 through photo-illumination. 
The dependence on U t was already discussed above; here 
the role of the other three parameters is investigated. 



1. Dependence on semiconductor width L a 

In Fig. QJJ we fixed the conductivity at its smallest 
value a s = 54/7000 = 7.714 • 10~ 3 and varied L s from 
10 to 150 in equal steps; resistance and capacitance then 
depend on L s like 1Z S oc L s and C s oc 1 / L s according to 
equation (|2ip in section Hi CI In physical units, the width 
of the semiconductor layer was changed from 0.28 mm to 
4.17 mm. The length of the gas gap was L — 36. Rather 
than the total applied voltage Ut, the current Jo = 1.32 • 
10~ 5 was fixed. For the value L s = 54, the corresponding 
parameter values are 1Z S = 7000, C s = 0.243 and Ut — 40 
as in Figs. [6] and [7] above. 

For small width L s the system shows a pure Hopf- 
instability (fc*=0), but with increasing L s , the most 



unstable mode k* suddenly becomes nonvanishing, i.e., 
the system undergoes a first order transition to spatio- 
temporal patterns. Furthermore, the growth rates Re A 
decrease and the oscillation frequencies Im A increase 
when L s increases while the Maxwell relaxation time 
t s = 1Z S C S of the semiconductor does not vary. 

For the smaller value 1Z S — 700, a similar behavior 
is observed: for the parameters of Fig. QJ there is still a 
Hopf bifurcation, but for increasing L s , the most unstable 
mode k* can become positive as well and a Turing- Hopf- 
instability occurs. 

2. Dependence on semiconductor resistance TZ 3 

Now the dependence on the conductivity of the semi- 
conductor is tested that can be varied by illumination 
by a factor of 10. Accordingly, in Fig. QTJ the resis- 
tance 1Z S varies in the interval between 700 and 7000, 
while the other parameters are chosen as in the previ- 
ous section: L = 36, L s = 54, C s = 0.243. The current 
Jo = 1.32 ■ 10~ 5 is fixed; this value corresponds to a volt- 
age of Ut = 40 for 1Z S = 7000. Therefore, the upper line 
in Fig. [IT] is the same as the one in Figs. OH [6] and [71 

For all values of the resistivity, the wave number k* 
where the growth rate Re A is maximal, stays positive 
(k* > 0), but it drops below zero for smaller 1Z S making 
the homogeneous stationary state stable. Furthermore, 
for the smallest investigated 1Z S , namely TZ S = 700, and 
fixed Jo, the voltage is U t = 16.25, but when U t = 23 
the growth rate is positive for the same k range and has 
maximum for k = as can be seen in Fig. 3] 

3. Dependence on length L of gas gap 

Fig- HH shows for the same parameter values that with 
increasing gas gap width L, the growth rate ReA(A;) in- 
creases and the oscillation frequency ImA(fc) decreases, 
while the most unstable mode stays nonvanishing: k* > 
0. The explored gap lengths L all correspond to a 
falling current voltage characteristics of the gas dis- 
charge, cf. section InFl 

V. NUMERICAL SOLUTIONS OF THE INITIAL 
VALUE PROBLEM 

A. Implementation and results 

The full dynamical problem was also solved numeri- 
cally as an initial value problem. This allows us to test 
the results of the stability analysis, to visualize the ac- 
tual dynamics and also to study the behavior beyond the 
range of linear stability analysis. Details of the numerical 
implementation are given in the appendix. 

We study the case of high resistivity 1Z S — 7000 that 
leads to spatial pattern formation as discussed above. 
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FIG. 11: The influence of the resistance lZ a of the semiconduc- 
tor layer on real and imaginary part of the dispersion relation 
A(fc) for equidistant 1Z S in the range between 1Z S = 700 and 
7000. The current is J = 1.32 • 10 -5 . All other parameters 
are as in Figs. H and [7] L = 36, L a = 54, C s = 0.243, and 
U t = 40 at lis = 7000. 
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FIG. 12: The influence of the width L of the gas gap on 
real and imaginary part of the dispersion relation X(k) for 
L = 30, 31, 32, 40. The current is J = 1.32 ■ 10 -5 . 

All other parameters are as in Figs. [5] [6] and [7] 1Z S = 7000, 
C s = 0.243, L s = 54, and U t = 40 at L = 36. 



multiple of the most unstable wave length 2ir/k* . After 
some tests with higher multiples showing essentially the 
same behavior, we used L x = 2 x 2ir/k*, i.e., the lateral 
extension is twice the expected wave length. The initial 
condition is 



p(x,z,0) = p (z) + C pi{z)e lk * x + 



(68) 



where k* is the wavenumber of the most unstable mode, 
Po(z) is the stationary solution, and p\{z) is the eigen- 
function for k — k* constructed in Sect. IIII CI The con- 
stant C is chosen such that the perturbation is small 
compared to Pq{z). Note that we need to specify the 
initial conditions for the ion density only, since electron 
density and field are determined instantaneously due to 
the adiabatic elimination of the fast electron dynamics. 

Figure [TBI shows about one period of oscillation within 
4 time steps for the pattern forming case (Ut = 46.4). 
For each instant of time, the rescaled electron density 
s = a I p and the ion density p are shown in the gas dis- 
charge region, and the electric field is shown both in the 
gas discharge and the semiconductor region, as will be 
discussed in more detail later. The upper row contains 
3D plots and the lower row contour plots. The figures 
show the characteristic electron and ion distribution of a 
glow discharge, but with a strong spatio-temporal mod- 
ulation. 

The temporal period predicted by linear perturbation 
theory is 528. This is agrees approximately with the nu- 
merical results. On the other hand, the destabilization 
of the homogeneous stationary state in Fig. [13] is already 
far developed and in the fully nonlinear regime. There- 
fore the results of the stability theory at this time give 
only an indication for the full behavior. In particular, the 
nonlinearity has created an onset to doubling the spatial 
period that is absent for small perturbations. 

For presenting the evolution in time, the spatial struc- 
ture has to be represented on a line rather than in 
the full (x, z)-plane. Obviously, the ion density on 
the semiconductor-gas-interface (x, L) is an appropriate 
quantity, since it characterizes the local intensity of the 
discharge glow in the transversal direction. Figures [14] 
and [l_5] show the complete temporal evolution in such 
a presentation. Fig. [TU presents data of a perturbation 
decaying towards the stationary homogeneous state for 
Ut — 23.7, while Fig. [TBI shows the growing destabiliza- 
tion of the homogeneous stationary state for Ut — 46.4; 
the late stage of this evolution was shown in Fig. [TBI 



B. Comparison of numerical and stability results 



Two values of the applied potential were investigated: 
Ut = 23.7 where the homogeneous stationary state is 
stable, and Ut = 46.4 where this state is Turing-Hopf- 
unstable. 

In the transversal direction, we use periodic bound- 
ary conditions. We choose the lateral extension L x as a 



When one wants to compare results of the numerical 
simulation and of the stability analysis, the evolution of 
different spatial modes has to be extracted from the sim- 
ulation. Appropriate quantities are the transversally av- 
eraged electric field on the gas-semiconductor interface 

£h(r) = ^- [ £ z {x,L,T)dx, (69) 

J-'x JO 
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and the spatial modulation of the field 

£ s (t) — max.£ z (x, L, t) — min£ z (x, L, r). (70) 

X X 

These quantities, or rather the logarithms of |6i(t) — 
Sq{L)\ and of £ s (r) are shown for the stabilizing case 
lit = 23.7 in Fig. [TBI and for the destabilizing case U t — 
46.4 in Fig. [TTJ 

In this logarithmic plot for the fields, the lines through 
the maxima are approximately straight, which means 



FIG. 13: Profiles and contour lines of electron and ion particle 
densities s = cr/fi and p in the discharge region, and electric 
field component E z in discharge and semiconductor region at 
time steps r = 7920, 8040, 8160, 8280 for U t = 46.4 and 
IZs = 7000. x and z coordinates are as in Fig. 1 and the text. 
For each time step, the data is represented as a 3D plot in 
the upper row and as a contour plot in the lower row. 
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FIG. 14: Evolution of ion density p(x, L, r) at the internal 
border z — L for Ut = 23.7 and 1Z S = 7000 as a function of 
the transversal coordinate x and time r. Note that time r 
increases towards the back. 




W IW 1 2D MD 330 3QD 33C' 4 DO <I90 

FIG. 15: The same as in the previous figure, but now for 
hit = 46.4; now the perturbations grow. Temporal snapshots 
in the full (x,z)-plane of the same numerical experiment are 
shown in Fig. 1131 




T 



FIG. 16: Temporal evolution of the transversally averaged 
electric field and of the spatial modulation of the field at the 
internal border z — L for Ut = 23.7 and 1Z S = 7000: (a) 
In \£h(i~) — £q(L)\ and (b) ui£ s (t) as a function of r. 




x 



FIG. 17: The same as in the previous figure, now for hit = 
46.4. 



that the growth is exponential. For the destabilizing case, 
the growth rate of the spatial mode £ s (r) is slightly larger 
than that of the homogeneous mode £^{t) which implies 
that the most unstable mode has a non-vanishing wave 
number k*: RcA(fc*) > ReA(0) in agreement with lin- 
ear stability analysis. Furthermore, at late stages when 
the dynamics is beyond the range of linear perturbation 
theory and becomes nonlinear, the growth of all modes 
accelerates. 
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FIG. 18: (Color online) Comparison of results of the PDE 
solutions (solid lines) and of the stability analysis (dashed 
line). Ion density p at the computational nodes between x — 
and x = L x /4 of the internal border z — L as a function of 
time for U t = 23.7 and TL S = 7000. 



Figs. [T5J and [lj|] show a quantitative comparison 
between stability analysis and computational results. 
Fig. [T5] shows the stabilizing case U t — 23.7. The sta- 
bility analysis predicts that k* = 0.050 is the most un- 
stable mode; it has the eigenvalue A(fc*) = —0.2807 • 
10~ 3 + 0.4320 • 10~ 2 i. Therefore, the period of the tem- 
poral oscillations is predicted as 27r/Im(A) = 1454, the 
characteristic decay time as l/Re(A) = 3563, and the 
characteristic wave length as 2n/k* — 126. This pre- 
dicted behavior is shown as the dashed line in the upper 
panel of Fig. [151 The solid lines show the numerical solu- 
tion, more precisely the time evolution of the ion density 
p(x, L, t) evaluated on the grid nodes in the range be- 
tween x = and x — L x /A on the gas-semiconductor 
interface z = L. Period and growth rate agree quantita- 
tively, therefore both simulations and stability analysis 
can be trusted. 

The predictions on the /c=0-mode are tested in the 
lower panel of Fig. [151 here the transversal extension 
of the simulation system was chosen so narrow that 
transversal modes had no space to develop: the width 
was taken as L x = 27r/(100£;*) where k* is the most un- 
stable wave length. In this case, only the k — 0-mode 
can grow, it has A(0) = -0.3547 ■ 10~ 3 + 0.7102 ■ 10~ 2 i. 
The plot again shows a very good agreement between 
stability analysis and simulation, now effectively for the 
one-dimensional case. 

Finally, in Fig. [T9l again the destabilizing state for Ut = 
46.4 is shown. The stability analysis predicts the most 
unstable wave number k* = 0.0267 and its eigenvalue 
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T 

FIG. 19: (Color online) The same as in the previous figure, 
now for Ut = 46.4. 



A(fc*) = 0.4615 • 10~ 3 + 0.1191 • lO-M. The two panels 
show again the predicted and the simulated oscillations 
in a laterally wide system allowing the formation of the 
fc*-mode, and in the narrow system that only has space 
for the fc=0-mode. Again, the agreement is convincing. 



VI. CONCLUSIONS AND DISCUSSIONS 

We have investigated the onset or decay of spatio- 
temporal patterns in a layered semiconductor-gas dis- 
charge system subject to a DC voltage. By means of 
linear stability analysis, we were able to derive complete 
phase transition diagrams, rather than only to investigate 
single points in parameter space in a simulation. 

We have used the simplest model possible: Only the 
drift motion of electron and ion densities is taken into 
account, and the electrons are adiabatically eliminated. 
The semiconductor is approximated as a linear Ohmic 
conductor, and nonlinear effects come in only through 
the space charges of the ions in the gas discharge gap and 
through the surface charges on the gas-semiconductor 
interface. (We remark that particle diffusion has a 
smoothening effect and is not expected to generate any 
new structures. However, effects like gas heating or non- 
linear semiconductor characteristics can create additional 
destabilization mechanisms.) 

Methods and results of the linear stability analysis of 
the homogeneous stationary state and of the full numeri- 
cal simulation of the initial value problem are presented. 
The choice ofparameters is guided by the experiment 
described in [7J; they are summarized in Sect. Ill El In 
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the experiment Q, the resistance of the semiconductor 
can be changed by a factor of 10 by photo-illumination 
without changing any other system parameter, and a full 
phase transition diagram is derived experimentally. It is 
seen that the system never relaxes to a spatially struc- 
tured time-independent state, but depending on the re- 
sistance, it either forms a homogeneous oscillating or a 
spatially structured oscillating state. 

Like in the experiments, the homogeneous stationary 
state is either stable, or it can be destabilized by a tem- 
poral (Hopf) or spatio-temporal (Turing-Hopf) mode; a 
purely spatial destabilization (Turing) is observed nei- 
ther in experiments nor in theory. The transition from 
stable to unstable behavior is in about the same param- 
eter regime in a diagram spanned by applied voltage Ut 
and inverse resistance of the semiconductor 1 /1Z S . How- 
ever, the parameter range where either Hopf- or Turing- 
Hopf-destabilization is predicted, does not agree with 
the experimental range where purely temporal or spatio- 
temporal patterns are observed. This is discussed in de- 
tail in Section TlV A 41 A possible explanation is that we 
might be comparing the wrong transitions. The theo- 
retical prediction of linear perturbation theory concerns 
a destabilization to weak running waves of long wave 
lengths. These waves resemble much more the "diffuse" 
moving bands reported only in the Ph.D. thesis 10], than 
the very nonlinear small dynamic filament structures de- 
scribed in [7}. This suggestion actually asks for fur- 
ther experimental investigations. It is interesting to note 
that the physical mechanism of these "diffuse" bands has 
nothing to do with particle diffusion, but only with the 
Laplacian nature of the created electric fields. For further 
predictions on the parameter dependence of the linear in- 
stabilities, we refer to Section HVBI 

Our numerical solutions of the initial value problem 
in Section [V] agree well with our linear stability analy- 
sis within its range of validity. First of all, this proves 
the correct implementation of both methods. Second, 
for larger amplitudes, new nonlinear spatial structures 
appear such as the spatial period doubling in Fig. 1131 
For TZ S = 7000 and Ut = 40, these oscillations actually 
have been seen to reach a limit cycle that corresponds to 
a standing wave. Of course, the full nonlinear behavior 
in three spatial dimensions should be investigated in the 
future, but the 2D results give a first indication for the 
behavior. 

We remark finally that only the simplest possible 
model with nonlinear space charge effects was investi- 
gated. Of course, the model can be extended by various 
additional mechanisms, but obviously the simple model 
already contains all relevant physics to predict the onset 
of pattern formation in the correct parameter regime. 
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APPENDIX A: NUMERICAL PROCEDURE 

Here we describe the numerical method used for solv- 
ing the initial value problem numerically. The compu- 
tation is based on a finite-difference technique to solve 
equations (EZJ) with boundary conditions (|2"g j) - (pjTj) 
and periodic boundary conditions in the transversal di- 
rection. 

The computational domain is a rectangular region 
[0, L x ] x [0, L z ] on a two-dimensional Cartesian coordinate 
system (x, z), which consists of two layers - gas and semi- 
conductor, see Figs. [Tl and |2"01 We use a uniform vertex- 
centered grid in the 'vertical' z-direction with nodes 

Zj =jAz, = j = 0,l,---,N 

and a uniform cell-centered grid with nodes 



A 

Ax = — , 
M 



1,2,- •• ,M 



on the 'horizontal' x-direction. The grid is spaced such 
that the internal interface between semiconductor and 
gas region lies exactly on the grid line. 

The densities s = a/fi and p and the electric potential 
(j) are evaluated on the nodes of the grid, while the x and z 
components of the electric field (£ x and £ z , respectively) 
are calculated on the surfaces of the computational cell 
(Fig. OP- 
TO obtain a finite-difference representation of equa- 
tions (|25|) and ([26]) . we first integrate them over the cell 
volume Xi-x/2 < x < x l+ i/ 2 , Zj-1/2 < z < z 3 +i/2- Let 
us consider in detail the equation for the ions ([2"o]) . After 
its integration, we come to 



dpj,i _ (£x P)j,i-l/2 ~ (£e P)j,i+l/2 

dr Ax 

(£ z p)j-l/2,i ~ (£z P)j+l/2,; 

Az 



+ fi,i .(Al) 



The subscripts i and j are related to x (transversal) and 
z (longitudinal) directions respectively and / stands for 
the source term of Eq. ([21?)) . 

The choice of Pj±x/2,i an d Pj,i±x/2 at the surfaces of the 
computational cell determines the concrete discretization 
method for the convective terms of Eq. (f2l)[) . We used the 
third-order upwind-biased scheme (see, e.g., [37| . p. 83), 
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FIG. 20: Computational domain and computational cell. 



which in z- and x-direction is given by 



l r 
6 



"z j+i/2,t i-Pj-u + 5 P'J^ + 2 Pj+i 

+£ zj+l/2,i ( 2 AjVi + 5 Pj+M - Pj+2,i) 
{£xP)j, i+ i/2 = g £ xj+l/2,i (~Pj,i-1- + 5 Pj,» + 2 Pi,i+l) 



(A2) 



Here, the electric field components are 



£ = max 



0, £ 



£- = 



0, £ 



and £ = — V0 is discretized as 



z j+l/2,i 



a; j,i+X/2 



Az 

i>j,i+l ~~ 4>j,i 



Ax 



(A3) 



For the numerical time integration, we used the ex- 
trapolated second order BDF2 method, see (37[, p. 204, 
[38j |. p. 197, whose variable step size version has the form 



(i+^r pm -i 



„m-2 



l + 2r 
(1 + r) 



l + 2r 



1 + 2r r 

Ar m {2F m - 1 - F m - 2 ) , m > 2, (A4) 



where the superscript m denotes the time r m with 
step size Ar m = r m — r r „_i and step size ratio r = 



Ar m /Ar m _i. Here F contains the discretized convective 
terms and a source term. Note that we have dropped spa- 
tial indices in Eq. (|A4|) . Since the two-step method needs 
p° and p 1 as starting values, the explicit Eulcr method 



Ar m F{r m - X ,p m - X ) 



(A5) 



is used for the first step m — 1 . Because of the explicit 
time integration, we are restricted by the standard CFL 
stability condition. 

The same space discretization technique and time in- 
tegration method are used also for the electron density 
equation (|25|) that contains no temporal derivative. Note 
that in this case the z-direction plays the role of 'time' 
in (TA~4l) and [KB. 



To obtain a finite-difference approximation for Pois- 
son's equation (|27|) , we use the traditional second order 
discretization: 



2$ 



(Ax) 2 
, 



P 



3,1 



(Az) 2 

gas-discharge layer, 
semiconductor layer. 



(A6) 



This equation is valid everywhere except at the gas 
semiconductor interface where one has to account for a 
finite surface charge as well as for a discontinuity of the 
dielectricity constant. On this interface, the discrete ver- 
sion of the 'jump' condition (|30[) is used instead of (|A6[) . 
The system of resulting difference equations is solved by a 
symmetrical successive over- relaxation method (SSOR), 
see [H, p. 343. 

The complete numerical procedure was organized as 
follows. For every new (m + l)th time step, first Pois- 
son's equation was solved using the known ion density p m 
and surface charge value q™ in the jump condition (|30| . 
determining the electric field components in the new time 
step. Second, the electron density in the new time step 
s m+i wag ca i cu i a tcd. This determined the source term 
in the continuity equation for ions. Third the ion density 
pm+i wag calculated, which finally determined the new 
value for the surface charge in ([317)) . 

The numerical convergence was checked by performing 
several calculations with different error tolerance param- 
eters for Poisson's equation, using refinement of the space 
grid, and different time stepping parameters. The num- 
ber of grid nodes used in the calculations was 52 x 361 
in the x and z directions, respectively, for the potential 
in the whole gas discharge - semiconductor region, and 
54 x 147 in the x and z directions for the particle densi- 
ties in the gas discharge region. When solving Poisson's 
equation, the iteration process is stopped when the rela- 
tive error is ||^ fe+1 ) - ^ \\/U (k+1) \\ < 5 • 1CT 7 . 
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